/* CFL_traits.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#include <src/miscellaneous/OS_basic.h>
#include <src/data/Operators.h>


#if !defined(_OS_CFL_)
#define _OS_CFL_
/*!
 Base behaviour for all CFL conditions
*/
template <class TimeOperator > 
struct CFL_traits {
	typedef typename TimeOperator						time_operator_type;
	typedef typename time_operator_type::argument_type	argument_type;
	typedef typename time_operator_type::time_type		time_type;
	typedef typename time_operator_type::result_type	result_type;
};


template <class TimeOperator > 
struct no_cfl : public CFL_traits<TimeOperator>
{
};
/*! \brief \f$ \frac{\Delta_t}{\Delta_x} < 1\f$

 classical transport equation CFL \f$ \frac{\Delta_t}{\Delta_x} < 1\f$ defined over a regular mesh on [0,1]
*/
template <class TimeOperator, std::size_t Dim = TimeOperator::Dim> 
struct Delta_t_Delta_x : public CFL_traits<typename TimeOperator> 
{
public : 
	OS_STATIC_CHECK(false); //to implement for your own set of Dimension
};
template <class TimeOperator>
struct Delta_t_Delta_x<TimeOperator,1> : public CFL_traits<typename TimeOperator> 
{
public : 
	Delta_t_Delta_x() : CFL_(1.){};
	time_type operator()(const time_operator_type& Op, argument_type S, const time_type & time, result_type & op_result ) const {
		return CFL_/S.size();
	}
	time_type CFL_;
};
template <class TimeOperator>
struct Delta_t_Delta_x<TimeOperator,2> : public CFL_traits<typename TimeOperator> 
{
public : 
	Delta_t_Delta_x() : CFL_(1.){};
	time_type operator()(const time_operator_type& Op, argument_type S, const time_type & time, result_type & op_result ) const {
		return CFL_/S.size(1);
	}
	time_type CFL_;
};
/*! brief \f$ \frac{\Delta_t}{\Delta_{x}^2} < 1\f$

 classical diffusive equation CFL \f$ \frac{\Delta_t}{\Delta_{x}^2} < 1\f$ defined over a regular mesh on [0,1]
*/
template <class T, std::size_t Dim = T::Dim> 
struct Delta_t_Delta_x_2 : public CFL_traits<T> 
{
public : 
	OS_STATIC_CHECK(false); //to implement for your own set of data
};
template <class T> 
struct Delta_t_Delta_x_2<T,1> : public CFL_traits<T> 
{
public : 
	Delta_t_Delta_x_2() : CFL_(1.){};
	time_type operator()(const time_operator_type& Op, argument_type S, const time_type & time, result_type & op_result ) const {
		return CFL_/S.size();
	}
	time_type CFL_;
};
template <class T> 
struct Delta_t_Delta_x_2<T,2> : public CFL_traits<T> 
{
public : 
	Delta_t_Delta_x_2() : CFL_(1.){};
	time_type operator()(const time_operator_type& Op, argument_type S, const time_type & time, result_type & op_result ) const {
		return (CFL_*CFL_/ ( S.size(1)*S.size(1) ));
	}
	time_type CFL_;
};

/* ! CFL condition for Positive defined operator
This CFL returns
\f[
    -2 \frac { \mathcal{O}(S) \cdot S} {|\mathcal{O}(S)|_{\ell^2}^2}
\f]

This CFL condition has been designed for the Crank_Nicolson class
working out positive defined operator ( \f$ \mathcal{O}(S) \cdot S < 0 \f$ ).
It ensure stability providing an uniform bound of the \f$ \ell^2\f$  norm of each iterate :
We recall that the fixed point algorithm for Crank Nicolson schemes writes

initialization : \f$ S^{0,n+1} : = S^{n}\f$. Loop for \f$k \ge 1\f$, compute the \f$k\f$-th iterate
\f[ \label{CN}
S^{k,n+1/2} := \frac {S^{n}+S^{k-1,n+1}}{2};S^{k,n+1} : = S^{n} + \tau^n \mathcal{O}(S^{k,n+1/2} )
\f]
Let \f$|S|_{\ell^2}^2 = \frac{1}{2} \sum |S|_i ^2\f$. We look for a \f$\ell_2\f$ stability of each iterate \f$S^{k,n}\f$ :
\f[
|S^{k,n+1}|_{\ell^2}^2 \le |S^{n}|_{\ell^2}^2
\f]

Taking the scalar product of the right-hand side of (\ref{CN}) with \f$ S^{k+1,n+1/2} = \frac {S^{n}+S^{k,n}}{2}\f$ leads us to
\f[
\frac {|S^{k,n+1}|_{\ell^2}^2 - |S^{n}|_{\ell^2}^2}{\tau^n} =  \mathcal{O}(S^{k,n+1/2} ) \cdot S^{k+1,n+1/2}
\f]


Equivalently we have stability if
\f[
\mathcal{O}( S^{k,n+1/2} ) \cdot S^{k+1,n+1/2} \le 0
\f]
We compute

\f[
S^{k+1,n+1/2} = S^{n} + \frac{\tau^n}{2} \mathcal{O}(S^{k,n+1/2} )
\f]

Thus we have stability if
\f[
 \mathcal{O}(S^{k,n+1/2} ) \cdot \left( S^{n} + \frac{\tau^n}{2} \mathcal{O}( S^{k,n+1/2} ) \right) \le 0
\f]

i.e.
\f[
\tau^n \le  -2 \frac { \mathcal{O}(S^{k,n+1/2}) \cdot S^{n} } {|\mathcal{O}(S^{k,n+1/2}|_{\ell^2}^2)}
\f]
*/
template <class TimeOperator> 
struct positive_define_cfl : public CFL_traits<typename TimeOperator> 
{
public : 
	time_type operator()(argument_type S, argument_type O_S ) const {
		OS_double norm22_op = jmmath::norm22(O_S);
		OS_double scalar_prod = jmmath::scalar_prod (S , O_S);
//! return \f$- \frac { Op(t^{1/2} ,S^{k,t+1/2) \cdot S^{t) ) {|Op(t^{1/2} ,S^{k,t+1/2)|_{\ell^2}^2}}\f$
		if (norm22_op > THRESHOLD)  {
			OS_DYNAMIC_CHECK (scalar_prod < 0., "positive_define_cfl : Operator is not positive defined. The cfl condition can not be computed in positive_define_cfl");						//works only for positive definite operator  
			return scalar_prod * ( -2. / (norm22_op)) ;	// ensure a conservation of the norm L2 : stability
		}
		return OS_MAX_REAL;
	}
};
template <class TimeTypeOperator> 
struct un_function : public operator_base<typename TimeTypeOperator::time_type, OS_double>
{
	un_function(){};
	OS_double operator()(argument_type time) const {
		return 1.;
	}
};

/*!  \brief CFL for bounded operators

    This CFL condition is for fixed point algorithm schemes using an operator 
    \f$\mathcal{O}\f$ operating over a collection of points \f$S^n=\{S_i^n\}_i\f$, the CFL is computed accordingly to
\f[
      \tau^n = \frac{C \epsilon}{ | \mathcal{O}(S^{n + 1/2})|_{\ell_2}^2 }
\f]
Default settings are
\f[
      C = 1.; \epsilon = 0.1
\f]

This CFL ensures stability providing an uniform bound of the (normalized) \f$\ell^2\f$ norm for operator \f$\mathcal{O}\f$ for which we know that the (normalized) scalar product
\f[
    \mathcal{O}(S) \cdot S = \sum_i \mathcal{O}(S)_i \cdot S_i \le C
\f]

    This case is useful for diffusive particle type operator. For instance, the multidimensional
    particle heat operator verifies \f$\mathcal{O}(S) \cdot S = 2D\f$, for \f$S = \{S_i \in \mathcal{R}^D\}_i\f$ .

    Under the previous assumptions, this C.F.L condition provides the bound
\f[
    |S^{k+1,n+1}|_{\ell_2}^2\le t^{n+1} (C+\epsilon)
\f]
   where we denoted the \f$k\f$-th iteration of the fixed point algorithm
\f[
    S^{k,n+1} = S^{n} + \tau^n \mathcal{O}(S^{k-1,n + 1/2})
\f]
and \f$S^{k-1,n + 1/2} = \frac{S^n+S^{k-1,n + 1}}{2}\f$.

    Proof : We compute
\f[
       \delta_t |S^{k,n}|_{\ell^2}^2 = S^{k,n + 1/2} \cdot
      \mathcal{O}(S^{k-1,n + 1/2})
\f]
Note that
\f[
    S^{k,n + 1/2} = \frac{S^{k,n+1} + S^n }{2} =  S^{k,n - 1/2} + \frac{\tau^n}{2} \mathcal{O}(S^{k-1,n + 1/2})
\f]

Thus we get
\f[
       \delta_t |S^{k,n}|_{\ell^2}^2 \le C + \frac{\tau^n}{2} | \mathcal{O}(S^{k-1,n + 1/2})|_{\ell_2}^2
\f]

that implies the bound \f$ \delta_t |S^{k,n}|_{\ell^2}^2 \le  C(1+\epsilon)\f$ under the assumtion
\f[
      \tau^n | \mathcal{O}(S^{k-1,n + 1/2})|_{\ell_2}^2 \le 2 \epsilon C
\f]

*/
template <class TimeOperator> 
struct controled_moment_cfl : public CFL_traits<TimeOperator>
{
public : 
	typedef typename operator_base<time_type, OS_double> Control_Function_type;
	controled_moment_cfl(OS_double C = 1) : C_(C)
	{};
	time_type operator()(argument_type OP_S) const {
		OS_double norm22_op = jmmath::moment2(OP_S);
		if (norm22_op > THRESHOLD)  {
			OS_double bound = C_*.1;
			return bound / norm22_op ;	// ensure a bound of the norm L2 : |S^{k+1,n+1/2)|_{\ell_2}^2 \le \sum f(t^{n+1/2})(1+.01)(t^{n+1/2}-t^{n-1/2}) \sim \int_0^{t^{n+1/2}}f(t)(1+.01)
		};
		return OS_MAX_REAL;
	}
private :
	OS_double C_;
};

/*!
This CFL condition ensure that each particles location after a a time step iteration
will remain in a small ball around its initial location : we want to ensure
|S^{k,n+1} - S^{t)|_{\ell^{\infty}} \le 0.2 * \inf_{i,j} |S_i^{t)-S_j^{t)|
and an iteration is 
S^{k,n+1}  = S^{t) + \Delta_t Op(t^{1/2} ,S^{k,n+1/2) ) ;
Thus we get 
\Delta_t \le 0.2 * \frac{\inf_{i,j} |S_i^{t)-S_j^{t)| }{|Op(t,S^{t))|_{\ell^{\infty}}}
*/
template <class TimeOperator> 
struct generic_particle_cfl : public CFL_traits<TimeOperator>
{
public : 
	time_type operator()(const time_operator_type& Op, argument_type S, const time_type & time, result_type & op_result = 0) const {
		OS_double norm_infty_op = jmmath::norm_infty(op_result);
		OS_double retuned_cfl (OS_MAX_REAL);
		if (norm_infty_op > THRESHOLD) {
			OS_double norm_infty_s = inf_distance(Op);
			retuned_cfl = 0.2* norm_infty_s/norm_infty_op;
		}
		return retuned_cfl;
	};
};


#endif
